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An  analysis  and  comparison  of  daily  and  yearly  solar  irradiation  from  the  satellite  CM  SAF  database  and  a 
set  of  301  stations  from  the  Spanish  SIAR  network  is  performed  using  data  of  2010  and  2011.  This  analysis 
is  completed  with  the  comparison  of  the  estimations  of  effective  irradiation  incident  on  three  different 
tilted  planes  (fixed,  two  axis  tracking,  and  north-south  horizontal  axis)  using  irradiation  from  these  two 
data  sources.  Finally,  a  new  map  of  yearly  values  of  irradiation  both  on  the  horizontal  plane  and  on  inclined 
planes  is  produced  mixing  both  sources  with  geostatistical  techniques  (kriging  with  external  drift,  KED). 

The  Mean  Absolute  Difference  (MAD)  between  CM  SAF  and  SIAR  is  approximately  4%  for  the  irradiation 
on  the  horizontal  plane  and  is  comprised  between  5%  and  6%  for  the  irradiation  incident  on  the  inclined 
planes.  The  MAD  between  KED  and  SIAR,  and  KED  and  CM  SAF  is  approximately  3%  for  the  irradiation  on 
the  horizontal  plane  and  is  comprised  between  3%  and  4%  for  the  irradiation  incident  on  the  inclined 
planes. 

The  methods  have  been  implemented  using  free  software,  available  as  supplementary  material,  and  the 
data  sources  are  freely  available  without  restrictions. 

©  2013  Elsevier  Ltd.  All  rights  reserved. 
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Nomenclature 

AEMET  Spanish  Meteorology  Agency 

/V 

pk  estimated  coefficients  of  the  deterministic  model  in 

kriging  with  external  drift 
BSRN  Baseline  Surface  Radiation  Network 
CM  SAF  Satellite  Application  Facility  on  Climate  Monitoring 

AG(0)  difference  between  GCMSAF(0 )  and  GSIAR( 0) 

ACgr(0)  difference  between  GKED( 0)  and  GCMSAF(0) 

AGred(O)  difference  between  GKED( 0)  and  GSIAR( 0) 
qCmsaf^ o)  yearly  global  irradiation  on  the  horizontal  plane  data 
obtained  by  estimations  from  CM  SAF 
Gsiar( 0)  yearly  global  irradiation  on  the  horizontal  plane  data 
obtained  by  on-ground  measurements  from  SIAR 

AGej rj<ED  difference  between  GKED( 0)  and  GCMSAF(0) 

A GsefARED  difference  between  GKED( 0)  and  GS,AR( 0) 

A Gef  difference  between  and  GSJAR 

e(s6)  interpolated  residual  in  kriging  with  external  drift 

Gked( 0)  yearly  global  irradiation  on  the  horizontal  plane 
estimated  with  kriging  with  external  drift 
Grfd( 0)  yearly  effective  global  irradiation  on  the  inclined 
plane  estimated  with  kriging  with  external  drift 
y(h)  semivariogram  function 

7(h)  estimator  of  the  semivariogram  function 


qcmsaf  yearly  effective  global  irradiation  incident  on  differ¬ 
ent  planes  estimated  from  data  from  CM  SAF 
GsJar  yearly  effective  global  irradiation  incident  on  differ¬ 
ent  planes  estimated  from  data  from  SIAR 
h  separation  vector  between  two  locations 

KED  kriging  with  external  drift 

Ai  kriging  weights  determined  by  the  spatial  depen¬ 

dence  structure  of  the  residual 
LUT  look-up  table 

m( Sq)  fitted  deterministic  part  of  the  random  spatial  field  at 

a  new  location 
MAB  Mean  Absolute  Bias 

MAD  Mean  Absolute  Difference 

MBD  Mean  Bias  Difference 

OK  ordinary  kriging 

qk( se)  auxiliary  predictors  obtained  from  the  fitted  values  of 

the  explanatory  variable  at  the  new  location  in  kri¬ 
ging  with  external  drift 
RMSD  biased  Root  Mean  Square  Difference 
RMSD*  unbiased  Root  Mean  Square  Difference 
RTM  radiative  transfer  model 

SIAR  Agroclimatic  Information  System  for  Irrigation 

SIS  shortwave  incoming  solar  radiation 

z  kriging  estimation  of  the  random  spatial  field 

Z(s)  random  spatial  field 


1.  Introduction 

Nowadays,  with  a  wide  range  of  applications  in  agriculture, 
climate  monitoring  and  renewable  energies,  research  in  solar 
irradiation  is  a  very  demanded  field. 

Solar  irradiation  can  be  evaluated  by  processing  images  from 
satellites  or  by  on-ground  measurements  with  pyranometers  in 
meteorological  stations.  The  high  cost  of  these  meteorological 
stations  and  the  requirement  of  specific  and  periodic  calibrations 
explain  the  low  density  of  the  existing  networks  in  many  countries, 
although  this  kind  of  measurements  is  reliable  to  elaborate  solar 
irradiation  maps  [1].  The  satellite  models  need  to  be  validated  and 
refined  with  high  quality  measurements,  which  are  provided  by  on¬ 
ground  stations  [2].  The  satellite  estimates  present  a  wide  spatial 
and  temporal  coverage,  but  their  spatial  resolution  is  in  the  range  of 
kilometres,  which  in  many  applications  may  not  be  sufficient  and 
this  can  be  improved  with  geostatistics  [3].  The  high  degree  of  site 
dependence  of  solar  irradiation  makes  geostatistics  suitable  to 
evaluate  the  spatial  distribution  of  solar  irradiation  and  to  build 
maps  with  pyranometers  measurements  [4,5]. 

Geostatistics  were  firstly  applied  in  the  study  and  estimation 
of  ore  resources  [6],  soil  properties  [7]  and  afterwards,  in  fields 
such  as  on-ground  water  analysis  [8]  and  solar  irradiation  maps 
with  kriging  techniques  [3,9].  Residual  and  ordinary  kriging  have 
been  applied  to  elaborate  solar  irradiation  maps  taking  into 
account  elevation  and  cloudiness  as  significant  variables  [10],  or 
topographic  shadow  cast  and  elevation  [11],  and  also  with 
artificial  neural  networks  (ANN)  with  temperature  and  precipita¬ 
tion  as  inputs  [12].  Kriging  with  external  drift  (KED)  has  been 
useful  to  develop  solar  irradiation  maps  using  multiple  linear 
regression  (MLR)  models  [13].  Comparing  solar  irradiation  maps 
obtained  with  different  techniques  and  inputs  is  necessary  to 
assess  the  divergence  of  the  estimates.  The  MESOR  project 
compared  EnMetSol,  Helioclim-2,  NASA  SSE  version  6,  Satel- 
Light  and  SOLEMI  databases  obtained  with  satellite  estimates 
and  ESRA,  PV  GIS  Europe,  and  Meteonorm  version  6.1  databases 


generated  from  geostatistical  models  and  meteorological  obser¬ 
vations  in  Europe  [1]. 

Recently,  the  Spanish  Agency  of  Meteorology  (AEMET)  has 
released  a  new  solar  irradiation  atlas  for  Spain  (the  former 
was  of  1984)  [14]  providing  monthly,  seasonal  and  annual 
average  of  global,  direct  and  diffuse  irradiation  on  the  hor¬ 
izontal  plane  with  a  resolution  of  3  km  using  monthly  data 
sets  from  1983  to  2005  of  CM  SAF.  Besides,  a  validation 
process  has  been  developed  comparing  CM  SAF  data  with 
uninterrupted  registers  from  2003  to  2005  of  29  meteorolo¬ 
gical  stations  from  the  National  Radiometric  Network  (RRN)  of 
AEMET.  On  the  other  hand,  for  direct  irradiation,  only  two 
ground  stations,  with  uninterrupted  data  from  1992  to  2005, 
were  selected.  The  Mean  Absolute  Deviation  (MAD)  obtained 
from  this  validation  process  for  global  monthly  average  is 
12.23  W/m2  (6.7%),  which  is  slightly  higher  than  the  CM  SAF 
target  of  10  W/m2.  It  is  important  to  underline  that  the  AEMET 
global  irradiation  atlas  is  restricted  to  the  horizontal  plane. 

This  paper  innovates  with  an  analysis  and  comparison  of  solar 
irradiation  from  the  CM  SAF  database  (Section  2.1)  and  a  large  set 
of  stations,  considering  301  meteorological  stations  (versus  the  29 
of  the  aforementioned  assessment  by  AEMET)  from  the  Spanish 
SIAR  network  (Section  2.2),  and  with  the  estimation  of  effective 
irradiation  incident  on  three  different  tracking  planes. 

Therefore,  the  contribution  of  this  paper  is  threefold: 

•  Analysis  and  comparison  of  daily  and  yearly  global  irradiation 
on  the  horizontal  plane  obtained  by  on-ground  measurements 
and  satellite  estimate  data. 

•  Analysis  and  comparison  of  yearly  global  irradiation  incident 
on  different  tilted  planes  (fixed,  two  axis  tracking  on  azimuth 
and  solar  elevation,  north-south  horizontal  axis)  estimated 
from  these  two  data  sources. 

•  Elaboration  of  a  new  map  of  yearly  values  of  irradiation  both 
on  the  horizontal  plane  and  on  inclined  planes  with  a  smooth 
combination  of  both  sources  using  geostatistical  techniques. 
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Fig.  1.  Hovmoller  plot  with  the  time  evolution  of  the  daily  horizontal  irradiation 
(Wh/m2)  as  published  by  CM  SAF,  averaged  along  10°W  to  5°E. 


The  analysis  comprises  daily  irradiation  data  of  2010  and 
2011.  The  global  irradiation  on  the  horizontal  plane  is  compared 
both  in  a  daily  and  a  yearly  basis,  while  the  effective  irradiation 
incident  on  different  planes  is  only  examined  in  a  yearly  basis.  In 
order  to  ease  the  discussion  of  results,  the  yearly  analysis  is 
carried  out  with  the  averages  of  2010  and  2011. 

To  enable  reproducible  research  [15],  the  methods  have  been 
implemented  using  free  software  (Section  6).  Both  the  source  code 
and  the  data  sources  are  freely  available  without  restrictions. 


2.  Radiation  data  sources 

2.2.  CM  SAF 

The  Satellite  Application  Facility  on  Climate  Monitoring  (CM 
SAF)  [16]  is  a  joint  venture  of  the  Royal  Netherlands  Meteorolo¬ 
gical  Institute,  the  Swedish  Meteorological  and  Hydrological 
Institute,  the  Royal  Meteorological  Institute  of  Belgium,  the 
Finnish  Meteorological  Institute,  the  Deutscher  Wetterdienst, 
Meteoswiss,  the  UK  MetOffice,  with  the  collaboration  of  the 
European  Organization  for  the  Exploitation  of  Meteorological 
Satellites  (EUMETSAT).  The  CM  SAF  was  funded  in  1992  to 
retrieve,  archive,  and  distribute  climate  data  to  be  used  for 
climate  monitoring  and  climate  analysis.  The  spatial  resolution 
of  the  different  products  ranges  from  15  km2  to  90  km2  [17]. 

The  CM  SAF  provides  two  categories  of  data:  operational 
products  and  climate  data.  The  operational  products  are  built  on 
data  that  is  validated  with  on-ground  stations  and  then  is 
provided  in  near  real  time  to  develop  variability  studies  in  diurnal 
and  seasonal  time  scales.  However,  climate  data  are  long-term 
data  series  to  assess  inter-annual  variability  [18]. 

In  this  study,  the  shortwave  incoming  solar  radiation  product 
(SIS)  is  selected  with  a  spatial  resolution  of  15  km2,  available  as 
daily  and  monthly  averages  (Figs.  1  and  2).  SIS  collates  shortwave 
radiation  (0.2-4  pm  wavelength  range)  reaching  an  horizontal 
unit  earth  surface  obtained  by  processing  information  from 
geostationary  satellites  (SEVIRI  sensor  on  board  of  the  METEOSAT 
Second  Generation  (MSG))  and  also  from  polar  satellites  (AVHRR 
sensor  on  NOAA  polar  satellites)  [17]  and  then  validated  with 
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Fig.  2.  Average  of  yearly  horizontal  irradiation  (kWh/m2)  on  the  horizontal  plane 
as  published  by  CM  SAF  during  2010  and  2011. 


Longitude 


Fig.  3.  Meteorological  stations  of  the  SIAR  network.  The  color  key  indicates  the 
altitude  (m).  Those  stations  whose  average  yearly  absolute  difference  from  the  CM 
SAF  values  is  higher  than  5%  are  displayed  with  red  points.  (For  interpretation  of 
the  references  to  color  in  this  figure  caption,  the  reader  is  referred  to  the  web 
version  of  this  article.) 


high-quality  on-ground  measurements  from  the  Baseline  Surface 
Radiation  Network  (BSRN).1 

In  this  paper,  SEVIRI  data  has  been  selected  following  the  CM 
SAF  recommendation  of  these  data  to  be  used  for  latitudes  south¬ 
ern  65°N  [19].  Validation  of  SEVIRI  SIS  data  with  4  BSRN  stations 
showed  that  more  than  90%  of  the  values  are  below  the  accuracy 
target  value  of  10W/m2  (plus  the  uncertainty  of  the  ground  based 
measurements).  Besides,  the  absence  of  a  trend  in  the  bias 
demonstrates  the  stability  and  homogeneity  of  the  product  [20]. 

The  method  for  retrieving  the  solar  surface  irradiance 
employed  by  CM  SAF  is  based  on  the  libRadtran  radiative  transfer 
model  (RTM)  [22]  in  combination  with  a  new  approach  of  several 
parameterizations  and  eigenvector  look-up  tables  (LUT).  A  LUT  is 
a  data  structure  with  discrete  pre-computed  RTM  results  for  a 
variety  of  atmospheric  and  surface  states.  Thus,  the  surface 
irradiance  (transmittance  multiplied  by  extraterrestrial  incoming 
solar  flux  density)  for  a  given  atmospheric  state  can  be  obtained 
by  interpolation,  through  the  LUTs,  for  each  satellite  pixel  and 
time.  Therefore,  with  a  LUT  approach  the  results  are  similar  to 
those  obtained  with  a  RTM  reducing  computation  costs  [23]. 


1  http:/ /www.bsrn.awi.de/en/home/. 
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Fig.  4.  Organization  of  the  analysis  procedure.  Ellipses  represent  point  data  sets  (values  from  the  meteorological  stations,  for  example)  and  rectangles  denote  raster  maps 
(values  from  CM  SAF,  for  example).  The  red  color  is  used  to  identify  the  original  sources,  green  for  comparison  results,  and  blue  for  transformation  results  (geostatistical 
interpolation  or  effective  irradiation),  (a)  Horizontal  irradiation  comparison,  (b)  Effective  irradiation  comparison.  (For  interpretation  of  the  references  to  color  in  this  figure 
caption,  the  reader  is  referred  to  the  web  version  of  this  article.) 


The  CM  SAF  method  still  can  be  improved  by  a  better  semi- 
empirical  adjustment  of  cloud  effects  and  by  improved  meteor¬ 
ological  information  about  aerosols  and  snow  cover  maps  [23]. 
In  fact,  one  main  goal  of  the  Continuous  Development  and  Operations 
Phase  of  the  CM  SAF  (2007-2012)  is  to  improve  all  data  sets  in  order 
to  develop  studies  of  inter-annual  variability  [17]. 

Fig.  1  displays  a  Hovmoller  plot  [21  ]  with  the  time  evolution  of 
CM  SAF  daily  irradiation  for  2010  and  2011  averaged  along  10°W 
to  5°E,  from  35.5°N  to  44°N.  Fig.  2  displays  the  average  of  annual 
global  irradiation  on  the  horizontal  plane  for  2010  and  2011. 

2.2.  SIAR 

Land-measured  daily  irradiation  is  collected  from  the  Agrocli- 
matic  Information  System  for  Irrigation  (SIAR)  [24]  a  free- 
download  database  operating  since  1999,  covering  the  majority 
of  the  irrigated  area  of  Spain  [25-29].  This  network  belongs  to  the 
Ministry  of  Agriculture,  Food  and  Environment  of  Spain,  as  a  tool 
to  predict  and  study  meteorological  variables  for  agriculture.  SIAR 
is  composed  by  12  regional  centers  and  a  national  center,  aiming 
to  centralise  and  depurate  measurements  from  the  361  stations  of 
the  network.  The  stations  include  SKYE-SP1110  (Campbell-Scien¬ 
tific)2  or  CMP6  (KIPP&ZONEN),3  first  class  pyranometers  according 
to  the  World  Meteorological  Organization  (WMO).4  The  absolute 
accuracy  is  within  ±  5%  and  is  typically  lower  than  ±  3%. 


2  ftp://ftp.campbellsci.com/pub/csl/outgoing/uk/manuals/spll  10.pdf. 

3  http://www.kippzonen.com/?product/1251/CMP  +  6.aspx. 

4  http://www.wmo.int/pages/index_es.html. 


The  calibration  of  the  pyranometers  is  performed  by  Tragsatec 
[24,30]  according  to  ISO  9847:1992  [31]  using  two  CMP6- 
KIPP&ZONEN  reference  pyranometers  [32]  on  a  yearly  basis. 
Irradiation  is  computed  on  a  half-hourly  basis  from  irradiance 
samples  recorded  each  10  s,  collated  through  a  CR10X  (Campbell 
Scientific)  datalogger  within  the  station  and  then  sent  to  the 
regional  and  national  centers  [24]. 

Data  has  been  filtered  under  two  assumptions:  average  annual 
irradiation  must  be  higher  than  1000  kWh/m2,  and  only  stations 
with  more  than  600  measurement  days  available  (out  of  a  total  of 
730,  2  years)  are  selected.  Besides,  some  stations  have  been 
omitted  due  to  difficulties  in  the  access  to  the  coordinates  of 
some  stations,  to  uncompleted  or  spurious  data  series,  or  to 
stations  out  of  the  area  of  study.  Eventually,  301  meteorological 
stations5  (Fig.  3)  and  their  daily  global  irradiation  measurements 
on  the  horizontal  plane  for  2010  and  2011  have  been  considered. 


3.  Methods 

3.2.  Statistics 

The  analysis  is  built  upon  the  next  structure  (Fig.  4): 

•  Analysis  and  comparison  of  daily  and  yearly  global  irradiation 
on  the  horizontal  plane  data  obtained  by  on-ground 


5  The  name  and  location  data  of  these  stations  are  available  at  http://solar. 
r-forge.r-project.org/data/SIAR.csv. 
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measurements,  GSIAR( 0),  and  satellite  data,  GCMSAF(0).  The  dif¬ 
ference  between  these  sources  is  a  matrix  of  values  examined 
in  Section  4.1. 

AG(0)  =  gcmsaf  (0) — Gsiar  (0)  (1 ) 


•  Analysis  and  comparison  of  yearly  global  irradiation  incident 
on  different  planes  (fixed,  two  axis,  north-south  horizontal 
axis)  estimated  from  these  two  data  sources,  GseIAR  and 
respectively.  The  difference  between  these  results  is  a  set  of 
three  matrices  examined  in  Section  4.1 


a  r  rCMSAF  rSIAR 
ALr ef  ~  uef  ~Uef 


/ - T  \  V2 

(ac2(0)) 

RMSDgo  -  v  J 

Gs,ar(  0) 

(8) 

„„„„  AC(0) 

MADG0  = 

Gs,ar(  0) 

(9) 

MBDC0-  AC(0) 

Gs,ar(  0) 

(10) 

3.2.  Effective  irradiation 

•  Elaboration  of  a  new  map  of  yearly  values  with  a  smooth 
combination  of  both  sources  using  geostatistical  techniques, 
both  for  the  horizontal  plane,  GK£D(0),  and  for  the  inclined 
planes,  GeFD.  Section  3.3  outlines  the  geostatistical  interpola¬ 
tion  technique  (kriging  with  external  drift,  KED)  used  to 
combine  the  information  from  SIAR  and  CM  SAF. 

The  difference  between  the  SIAR  stations  and  the  results  of  the 
interpolation  are 

agk£d(°)  =  Gked(0)-Gsiar(0)  (3) 


A  rSlAR  _  rKED 
ALre/,/CED  ~  u ef 


rSIAR 

e/ 


The  difference  between  the  CM  SAF  maps  and  the  results  of  the 
interpolation  are 


A  rCMSAF 
ALrKED 


(0)  =  gked 


(0)-Gcmsaf(0) 


A  rCMSAF  _  rKED  rCMSAF 
ALje/,/C£D  —  Uef  ~Uef 


These  differences  are  summarized  using  several  statistics:  the 
unbiased  and  biased  Root  Mean  Square  Difference  (RMSD*  and  RMSD, 
respectively),  the  Mean  Bias  Difference  (MBD)  and  the  Mean  Absolute 
Difference  (MAD)  (Tables  3-5).  It  must  be  noted  that  RMSD2  = 
RMSD*2 -h  MBD2  and  that  MAD  <  RMSD  <  n1/2  -MAD.  The  reader  is 
referred  to  Ref.  [33]  for  the  details  on  the  convenience  of  these 
statistics. 

These  statistics  are  normalized  by  the  average  SIAR  values 
when  comparing  CM  SAF  or  KED  with  SIAR  (Eqs.  (l)-(4))  and  by 
the  average  CM  SAF  values  when  comparing  KED  with  CM  SAF 
(Eqs.  (5)  and  (6)).  For  example,  the  RMSD*,  RMSD,  MBD  and  MAD 
corresponding  to  Eq.  (1)  are 


RMSD*g  o  = 


qCMSAF ^  _  qCMSAF 


(0))  -  (csmr(0)-Gsmr(0))]  2 ) 


1/2 


qSiar  (0) 


(7) 


Three  different  tracking  methods  have  been  considered 

•  Fixed  plane,  oriented  towards  the  south  and  with  optimum 
inclination  angle  (latitude  minus  10°). 

•  North-south  horizontal  axis  tracker:  the  axis  of  rotation  is 
horizontal  with  respect  to  the  ground  and  is  on  a  north-south 
line.  Panels  are  mounted  horizontally  upon  the  tube  which  will 
rotate  on  its  axis  to  track  the  apparent  motion  of  the  sun 
through  the  day. 

•  Two-axis  tracker:  both  azimuth  and  altitude  are  constantly 
changing  to  track  the  sun. 

Detailed  description  of  these  methods  can  be  found  in  [34]. 
Table  1  summarizes  the  calculation  procedure  from  global  daily 
irradiation  on  the  horizontal  plane  to  effective  global  irradiation 
incident  on  an  inclined  plane.  It  is  important  to  highlight  that  this 
calculation  procedure  does  not  include  shadow  losses. 

The  first  step  of  the  procedure  (once  sun  and  trackers  geome¬ 
try  equations  have  been  computed)  is  to  decompose  the  daily 
global  irradiation  on  the  horizontal  plane  in  two  components, 
direct  and  diffuse  irradiation.  Diffuse  fraction,  the  ratio  of  diffuse 
to  global  irradiation,  is  estimated  from  the  clearness  index  with 
the  equations  proposed  in  [35]. 

The  second  step  is  to  build  irradiance  profiles  from  daily 
irradiation  values.  The  ratio  of  the  diffuse  irradiance  to  diffuse 
irradiation  is  assumed  to  be  equivalent  to  the  ratio  of  extra¬ 
terrestrial  irradiance  to  extraterrestial  irradiation.  The  ratio  of 
global  irradiance  to  daily  global  irradiation  is  estimated  with 
the  equations  proposed  in  [35].  It  must  be  noted  that,  because  of 
the  frequent  low  variability  of  solar  irradiance,  this  step  assumes 
that  the  average  value  of  irradiance  during  a  short  time  interval 
(for  example,  an  hour)  coincides  numerically  with  the  irradiation 
during  that  interval.  Under  this  assumption  the  profile  of  irra¬ 
diance  incident  on  a  surface  estimated  in  the  next  step  can  be 
aggregated  to  produce  daily  irradiation. 

The  third  step  computes  direct  and  diffuse  irradiance  incident 
on  the  inclined  plane  considering  purely  geometrical  criteria. 
Direct  irradiance  is  estimated  with  the  solar  zenith  angle  and 
the  angle  of  incidence  on  the  generator.  Diffuse  irradiance  is 
calculated  with  the  anisotropic  model  proposed  in  [36].  This 


Table  1 

Calculation  procedure  for  the  estimation  of  effective  irradiation  incident  on  a  PV  generator  from  daily  global  horizontal  irradiation  data. 


Step 


Method 


Sun  and  trackers  geometry 
Decomposition  of  daily  global  horizontal 
irradiation 

Estimation  of  irradiance 

Estimation  of  irradiance  on  inclined  surface 

Albedo  irradiance 

Effects  of  dirt  and  angle  of  incidence 


Set  of  equations  as  provided  in  [34] 

Correlation  between  diffuse  fraction  of  horizontal  irradiation  and  clearness  index  [35] 

Ratio  of  global  irradiance  to  daily  global  irradiation  [35] 

The  direct  irradiance  is  calculated  with  geometrical  equations.  The  estimation  of  the  diffuse  component  makes  use  of  the 
anisotropic  model  [36] 

Isotropic  diffuse  irradiance  with  reflection  factor  equal  to  0.2. 

Equations  proposed  in  [37].  A  low  constant  dirtiness  degree  has  been  supposed  (2%) 
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model  divides  the  diffuse  irradiance  in  circumsolar  (near  the  sun 
region)  and  isotropic,  using  an  anisotropy  index  to  estimate  the 
ratio  between  them.  The  albedo  is  assumed  to  be  isotropic  and  is 
estimated  from  the  global  irradiance  with  a  reflection  factor 
of  0.2. 

The  last  step  estimates  the  effective  irradiance  incident  on  a 
generator  subtracting  dust  and  angle  of  incidence  losses  from  the 
incident  irradiance  with  the  model  proposed  in  [37]. 

3.3.  Geostatistical  interpolation 

Geostatistics  deals  with  the  analysis  of  random  fields  Z(s), 
where  s  =  (x,y)  is  a  location  and  x,  y  its  geographical  coordinates. 
Measurements  of  the  random  field  Z  is  commonly  only  available 
at  a  limited  set  of  locations  (in  our  case,  the  meteorological 
stations  of  the  SIAR  network).  In  order  to  predict  the  value  of  Z  at 
locations  without  observations,  the  geostatistical  analysis 
involves  estimation  and  modelling  of  spatial  correlation  under 
simplifying  assumptions  of  stationarity  [38]. 

Assuming  that  the  samples  are  representative,  non- 
preferential  and  consistent,  values  of  the  field  at  a  new  location 
can  be  derived  using  a  spatial  prediction  model.  This  geostatistical 
interpolation  procedure  is  generally  known  as  kriging  [39].  A 
standard  version  of  kriging  is  called  ordinary  kriging  (OK).  Here 
the  predictions  are  based  on  the  model: 

Z(S)  =  jU  +  £(S)  (11) 

where  p  is  the  constant  stationary  function  (global  mean)  and  e(s) 
is  the  spatially  correlated  stochastic  part  of  variation. 

The  assumption  of  constant  mean  is  hardly  acceptable  for 
the  estimation  of  irradiation  over  a  large  area.  Ordinary  kriging 
was  initially  tried  within  this  study,  generating  inaccurate  esti¬ 
mations  due  to  the  long  distances  among  some  stations.  In 
mountainous-heterogeneous  regions  such  as  Galicia  (north  of 
Spain),  this  inaccuracy  was  more  significant  than  in  flat- 
homogeneous  regions,  such  as  Castilla-La-Mancha  (center  of 
Spain). 

This  model  can  be  improved  including  additional  information 
from  an  exhaustively-sampled  explanatory  variable.  If  the  expla¬ 
natory  variable  is  significantly  correlated  with  the  field  Z(s), 
predictions  at  a  new  location,  se,  can  be  obtained  modelling  the 
deterministic  and  stochastic  components  separately 

z(s0)  =  m(s0)  +  e(s0)  (12) 

where  m(se)  is  the  value  of  the  fitted  deterministic  part  at  the 
new  location,  c(Sq)  is  the  interpolated  residual.  These  two  com¬ 
ponents  can  be  derived  with 

Z(Sg)  =  Y)  PkQki s»)+  J2  (13) 

k= 0  1=1 

/V 

where  are  the  estimated  coefficients  of  the  deterministic 
model,  qk( s0)  are  the  auxiliary  predictors  obtained  from  the  fitted 
values  of  the  explanatory  variable  at  the  new  location,  2/  are  the 
kriging  weights  determined  by  the  spatial  dependence  structure 
of  the  residual,  and  e(s,)  are  the  residual  at  location  s*. 

This  improved  model  (Eq.  (13))  is  known  as  kriging  with 
external  drift  (KED)  or  regression  kriging  [39].  In  this  paper,  the 
explanatory  variable  is  the  irradiation  on  the  horizontal  plane 
estimated  by  CM  SAF,  gcmsaf( 0),  both  for  the  irradiation  on  the 
horizontal  plane  and  for  the  irradiation  incident  on  inclined 
planes.  Therefore,  the  KED  method  is  fed  with  three  sources  of 
information  to  produce  new  maps: 

•  Yearly  irradiation  measurements  on  the  horizontal  plane  from 

SIAR  stations  or  estimations  of  yearly  irradiation  on  the 

inclined  plane  based  on  the  measurements  from  SIAR. 


Table  2 

Parameters  of  the  variograms  fitted  to  the  SIAR  data  for  different  planes  using  CM 
SAF  irradiation  as  explanatory  variable. 


Irradiation  plane 

Model 

Nugget 

Sill 

Range 

C(  0) 

Pure  nugget 

4609.11 

— 

— 

Fixed 

Pure  nugget 

7275.30 

- 

- 

N-S  horizontal 

Spherical 

13,138.08 

6768.78 

458.45 

Two  axis 

Spherical 

19,831.10 

9336.59 

478.47 

Fig.  5.  Semivariances  and  variograms  fitted  to  the  SIAR  data  for  different  planes 
using  CM  SAF  irradiation  as  explanatory  variable. 

•  Estimations  of  yearly  irradiation  on  the  horizontal  plane  from 
CM  SAF  as  explanatory  variable. 

•  A  semivariogram  function  to  model  the  spatial  dependence 
structure  of  the  residuals. 

The  semivariogram  is  a  function  defined  as  [40,41  ] 

7(h)  =  ly.(£(s)-e(s + h))2  (14) 

where  h  is  the  separation  vector  between  two  locations,  h  =  Sj-Sj. 
This  equation  is  defined  under  the  assumption  that  the  variance 
of  e  is  constant  and  that  spatial  correlation  of  e  does  not  depend 
on  location  s  but  only  on  separation  distance  h.  The  estimator  of 
the  variogram,  called  the  sample  semivariogram  is 

y(h)  =  ,15) 

h  N, 

with  Nh  =  {(s,-,Sj) :  s—Sj  =  h},  the  set  of  all  pairs  of  locations 
separated  by  vector  h.  It  is  common  to  assume  that  the  variogram 
is  isotropic  and,  consequently,  that  the  correlation  at  two  loca¬ 
tions  depends  only  on  the  distance  between  them  and  not  on  the 
direction  between  them. 

The  sample  variogram  gives  estimates  only  at  observed  spatial 
lags.  Therefore,  it  is  not  enough  for  prediction  at  new  locations.  A 
common  solution  is  to  infer  a  parametric  variogram  model  from 
the  data  fitting  a  model  to  the  sample  variogram.  Some  well- 
known  parametric  variogram  functions  are  the  exponential, 
gaussian  or  spherical  models.  The  parameters  of  the  model  to 
be  determined  are  the  sill,  the  range  and  the  nugget  [39].  Table  2 
displays  the  parameters  of  the  variograms  fitted  to  the  SIAR  data 
for  different  planes  using  CM  SAF  irradiation  as  explanatory 
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a 


b 


GO 


Fixed 


[-0.12,-0.032) 

[-0.032,0.0038) 

[0.0038,0.031) 

[0.031,0.06) 

[0.06,0.095) 

[0.095,0.15) 

[0.15,0.27) 

[0.27,0.41] 


Fig.  6.  Relative  differences  of  the  yearly  values  of  horizontal  (Eq.  (1))  and  effective  irradiation  (Eq.  (2))  between  CM  SAF  and  SIAR  for  the  whole  set  of  SIAR  stations, 
(a)  Dotplot  and  (b)  Map. 


Table  3 

Statistics  of  the  yearly  irradiation  values  from  CM  SAF  and  SIAR.  The  RMSD,  RMSD*,  MBD  and  MAD  statistics  are  calculated  with  adimensionalized  differences  using  GsyIAR( 0) 
or  GseIAR  as  normalization  factors. 


Irradiation  plane 

Ocmsaf  (l<Wh/m2) 

& siar  (l<Wh/m2) 

MBD  (%) 

RMSD *  (%) 

RMSD  (%) 

MAD  (%) 

GO 

92.50 

102.09 

3.41 

4.44 

5.60 

4.19 

Fixed 

88.16 

112.09 

3.59 

5.21 

6.33 

4.69 

N-S  horiz 

146.21 

170.28 

4.24 

5.93 

7.30 

5.57 

Two-axis 

155.99 

195.63 

4.33 

6.36 

7.69 

5.84 
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Time 


Fig.  7.  Time  evolution  of  the  relative  differences  between  the  daily  global 
irradiation  on  the  horizontal  plane  from  SIAR  and  CM  SAF.  The  red  line  represents 
the  median  and  the  blue  lines  represent  the  5%  and  95%  quantiles. 
(For  interpretation  of  the  references  to  color  in  this  figure  caption,  the  reader  is 
referred  to  the  web  version  of  this  article.) 


variable.  Fig.  5  shows  the  semivariances  and  the  fitted  variogram 
models. 

The  nugget  effect,  associated  to  micro-variability  and  mea¬ 
surement  error,  models  the  discontinuity  of  the  variogram  at  the 
origin.  When  the  nugget  effect  is  present,  the  kriging  method  is 
not  an  exact  interpolator  (it  does  not  preserve  the  original 
observations).  It  must  be  highlighted  that  the  variograms  corre¬ 
sponding  to  irradiation  on  the  horizontal  plane  and  on  a  fixed 
plane  are  the  pure  nugget  model,  that  is,  the  residuals  show  no 
spatial  auto-correlation. 


4.  Discussion  of  the  results 


4.2.  Comparison  between  SIAR  and  CM  SAF 


The  comparison  of  GCMSAF(0)  and  GSIAR( 0)  must  be  performed 
taking  into  account  that  the  SIAR  pyranometers  present  a  toler¬ 
ance  of  5%  (Section  2.2).  In  Fig.  6a  71%  of  the  locations  are  inside 
the  range  of  this  pyranometer  uncertainty.  Outside  this  5%  band, 
96.5%  of  the  stations  SIAR  provide  lower  global  irradiation  values 
than  CM  SAF. 


The  relative  difference  increases  when  a  tracking  system  is 
considered,  |  AG(0)|/GSMr(0)  <  |  AGef\/GeftSiAR  (Fig.  6).  Besides, 
Table  3  shows  that  both  the  standard  deviation  of  the  irradiation 
values  {a siar  and  &cmsaf)  and  the  statistics  of  the  differences 
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Fig.  8.  Statistics  of  the  daily  global  irradiation  on  the  horizontal  plane  from  SIAR  and  CM  SAF.  (a)  Absolute  value  of  the  MBD.  (b)  MAD.  (c)  RMSD  and  (d)  RMSDc. 
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Fig.  9.  Global  solar  irradiation  estimated  with  KED  using  CM  SAF  as  external  drift,  (a)  GO.  (b)  Fixed,  (c)  NS  Horiz.  (d)  Two. 
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Fig.  10.  Histograms  of  the  normalized  differences  between  the  effective  irradia¬ 
tion  incident  on  a  fixed  plane,  a  north-south  horizontal  axis  tracker,  and  a  two- 
axis  tracker. 


i  — i - • - r - - - 1 

-10  -50  5  10 

Fig.  11.  Relative  differences  (%)  between  the  horizontal  irradiation  estimated  with 
kriging  of  the  values  at  the  SIAR  stations  using  the  Gy(0)  from  CM  SAF  as  external 
drift,  and  the  horizontal  irradiation  estimated  from  CM  SAF.  Positive  values  mean 
that  the  estimation  with  kriging  is  higher  than  with  CM  SAF. 


Table  4 

Statistics  of  the  yearly  horizontal  irradiation  values  from  KED  and  SIAR.  The  RMSD, 
RMSD*,  MBD  and  MAD  statistics  are  calculated  with  adimensionalized  differences 
using  GsyIAR( 0)  or  GseIAR  as  normalization  factors. 


Irradiation 

&  KED 

G  SIAR 

MBD 

RMSD * 

RMSD 

MAD 

plane 

(l<Wh/m2)  (l<Wh/m2) 

(%) 

(%) 

(%) 

(%) 

GO 

72.67 

102.09 

0.00 

4.28 

4.28 

3.02 

Fixed 

63.01 

112.09 

0.00 

5.23 

5.23 

3.66 

N-S  horiz 

122.67 

170.28 

0.00 

4.60 

4.60 

3.32 

Two-axis 

131.11 

195.63 

0.00 

5.01 

5.01 

3.63 

( RMSD  and  MAD )  increase  with  the  application  of  the  formulas  to 
account  for  tilted  surfaces.  This  observation  is  consistent  with 
[42,43],  where  the  variability  of  the  effective  irradiation  incident 
on  tracking  planes  was  reported  to  be  higher  than  the  variability 
of  irradiation  on  the  horizontal  plane. 

No  significant  latitudinal  behavior  is  appreciated  in  any  of  the 
cases  of  Fig.  6,  although  as  per  Figs.  1  and  2,  solar  irradiation  is 
clearly  latitudinally  dependent. 

In  Fig.  7,  AG(0)  presents  a  seasonal  periodicity  of  the  5%  and  95% 
quantiles,  with  a  wider  range  for  winter  and  more  confined  in 
summer.  In  this  figure  SIAR  presents  a  set  of  samples  in  which 
Gsiar( 0)  is  significantly  lower  than  gcmsaf( 0).  It  may  be  explained 
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Fig.  12.  Relative  differences  of  the  yearly  values  of  horizontal  (Eq.  (3))  and  effective  irradiation  (Eq.  (4))  between  KED  and  SIAR  for  the  whole  set  of  SIAR  stations, 
(a)  Dotplot.  (b)  Map. 
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Fig.  13.  Relative  differences  (%)  between  the  effective  irradiation  incident  on  a 
fixed  plane  estimated  with  kriging  of  the  values  at  the  SIAR  stations  using  the 
Gy(0)  from  CM  SAF  as  external  drift,  and  the  effective  irradiation  estimated  with 
the  CM  SAF  raster.  Positive  values  mean  that  the  estimation  with  kriging  is  higher 
than  with  CM  SAF. 


Fig.  15.  Relative  differences  (%)  between  the  effective  irradiation  incident  on  a 
two-axis  tracker  estimated  with  kriging  of  the  values  at  the  SIAR  stations  using  the 
Gy(0)  from  CM  SAF  as  external  drift,  and  the  effective  irradiation  estimated  with 
the  CM  SAF  raster.  Positive  values  mean  that  the  estimation  with  kriging  is  higher 
than  with  CM  SAF. 
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Fig.  14.  Relative  differences  (%)  between  the  effective  irradiation  incident  on  a 
north-south  horizontal  axis  tracker  estimated  with  kriging  of  the  values  at  the 
SIAR  stations  using  the  Gy(0)  from  CM  SAF  as  external  drift,  and  the  effective 
irradiation  estimated  with  the  CM  SAF  raster.  Positive  values  mean  that  the 
estimation  with  kriging  is  higher  than  with  CM  SAF. 


due  to  local  events  not  registered  by  the  satellite  resolution,  or  to 
failures  in  the  on-ground  registers,  which  were  not  detected  when 
filtering  spurious  data.  It  is  important  to  highlight  that  these 
extreme  events  are  smoothed  with  the  averages  of  annual  sums. 

In  Fig.  8a  and  b,  the  statistics  MBD  and  MAD  are  lower  than  5% 
in  most  of  the  stations,  although  a  set  of  outliers  is  appreciated  in 
the  Valencia  region  (middle  east  of  Spain).  In  Fig.  8c  and  d,  the 
RMSD  and  RMSD*  are  generally  lower  than  7%.  In  a  set  of  stations  in 
the  north  of  Spain  in  which  the  MBD  were  lower  than  6%,  the  RMSD 
and  RMSD *  are  significantly  higher,  which  may  be  explained  due  to 
the  strong  meteorological  variability  existing  in  the  Ebro  valley. 
In  the  middle  Ebro  valley  there  are  marked  thermal  contrasts, 
with  possible  generation  of  orographic  fog,  typical  in  valleys. 


This  variability  generates  a  much  more  variable  distribution  of 
error  magnitudes  which  produces  higher  levels  of  RMSD  [33]. 

Both  Figs.  6  and  8  compare  irradiation  on  the  horizontal  plane 
from  SIAR  and  CM  SAF.  However,  there  are  remarkable  differences 
between  them.  For  example,  there  are  some  stations  in  the  north 
of  Spain  clearly  visible  in  Fig.  6  (important  difference  between  CM 
SAF  and  SIAR)  but  they  are  invisible  in  Fig.  8.  To  explain  this 
apparently  contradictory  behavior  is  important  to  note  that  some 
stations  include  missing  values  in  their  data  sets.  Fig.  8  compares 
daily  values  with  a  collection  of  statistics  computed  without 
those  missing  values.  However,  Fig.  6  compares  yearly  values 
with  missing  values  contributing  as  zeros.  Therefore,  those  sta¬ 
tions  with  a  higher  proportion  of  missing  values  will  provide 
lower  annual  irradiation  values,  although  their  daily  statistics 
could  be  assimilable  to  a  station  without  missing  values. 

4.2.  Comparison  between  I<ED,  SIAR  and  CM  SAF 

The  KED  technique  does  not  perform  as  an  interpolation 
function  when  the  nugget  effect  is  not  null,  which  occurs  in  this 
study  as  shown  in  Table  2.  This  fact  indicates  that  there  is  an 
intrinsic  variability  independent  from  the  distance  between  sta¬ 
tions.  In  this  case,  the  KED  behaves  as  a  smoothing  function  of  the 
SIAR  values,  with  the  external  drift  of  CM  SAF,  generating  a 
solution  that  differs  both  from  SIAR  and  CM  SAF  in  the  positions 
of  the  meteorological  stations. 

Maps  of  global  irradiation  obtained  with  KED  using  CM  SAF  as 
external  drift  are  shown  in  Fig.  9  on  the  horizontal  surface,  fixed 
tilted  plane,  tracking  system  with  north-south  axis  and  two-axis 
tracking  system.  Differences  between  one  axis  tracking  and  fixed 
plane  range  from  0.2%  to  36%,  and  between  two  axis  tracking  and 
fixed  plane  range  from  11%  to  55%  (Fig.  10).  These  differences  are 
more  significant  in  southern  Huesca  (42°N,  0°W),  Zamora  (42°N, 
6°W),  the  peninsular  center  (38°N  to  41  °N,  1°W  to  6°W)  and 
Almeria  (37°N,  2°W),  and  lower  along  the  Cantabric  coast  (43°N, 
1°W  to  9°W),  due  to  the  reduced  influence  of  direct  irradiation. 
Differences  in  irradiation  estimated  with  KED  between  the  two 
tracking  systems  range  from  1 1%  to  14%,  with  higher  values  in  the 
Ebro  valley  (40°N  to  42°N,  5°W  to  2°E)  and  along  the  Mediterra¬ 
nean  coast  and  lower  values  in  Jaen  (38°N,  4°W). 
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Fig.  16.  Violin  plot  of  the  relative  differences  (%)  between  irradation  estimated  with  kriging  of  the  values  at  the  SIAR  stations  using  CM  SAF  as  external  drift,  and  the 
irradiation  estimated  with  the  CM  SAF  raster.  Positive  values  mean  that  the  estimation  with  kriging  is  higher  than  with  CM  SAF.  Each  latitude  interval  include  between  41 
and  43  stations  with  an  overlap  of  10%. 


Table  5 

Statistics  of  the  yearly  horizontal  irradiation  values  from  KED  and  CM  SAF  (Figs.  1 1 
and  13-15).  The  RMSD,  RMSD*,  MBD  and  MAD  statistics  are  calculated  with 
adimensionalized  differences  using  GyMSAF( 0)  or  G™SAF  as  normalization  factors. 


G  KED 

(l<Wh/m2) 

G  CMSAF 

(l<Wh/m2) 

MBD 

(%) 

RMSD* 

(%) 

RMSD 

(%) 

MAD 

(%) 

GO 

112.24 

142.86 

-2.52 

1.84 

3.12 

2.86 

Fixed 

93.06 

148.43 

-2.31 

2.93 

3.73 

3.42 

N-S  horiz 

215.48 

239.73 

-3.48 

2.62 

4.35 

3.60 

Two-axis 

242.62 

265.95 

-3.67 

2.75 

4.59 

3.82 

Sked(  0) 


Figs.  11,  6b  and  12b  reveal  that  AG 
A-GSJfAKED\/Gef,siAR  are  slightly  lower  than  |AC(0) 


/Gsiar( 0)  and 
/Gsmf(O)  and 

A Gef  /GefjSIAR,  respectively,  and  correspondingly  for  the  CM  SAF 
values.  Once  again,  the  RMSD  and  MAD  values  are  very  similar  for 
all  trackers  and  higher  than  those  corresponding  to  the  irradiation 
on  the  horizontal  plane  (Table  4). 

In  Fig.  16,6  higher  latitudes  present  higher  dispersion  of  the 
differences  than  lower  latitudes,  although  values  remain  in  a  4% 
band.  Specially,  from  40°N  to  north,  just  when  average  elevation 
increases  (Fig.  3),  dispersion  values  are  higher.  As  already  men¬ 
tioned,  CM  SAF  shows  a  more  inaccurate  behavior  when  clouds  or 
snow  can  appear.  This  fact  can  widen  the  range  of  differences  for 
mountainous  areas.  In  Fig.  2  mountainous  areas  act  as  modulators 
of  irradiation  [14].  In  Fig.  11  (irradiation  on  the  horizontal  plane) 
SIAR  only  presents  higher  irradiation  than  CM  SAF  in  the  very 
north  of  Spain.  Nevertheless,  the  variability  of  the  previous  map  is 
in  the  range  of  5%,  which  stands  within  the  uncertainty  band. 


6  This  figure  displays  the  data  distribution  with  a  violin  plot,  a  combination  of 
a  boxplot  and  a  kernel  density  plot.  Therefore,  this  graphical  tool  shows  the  lower 
quartile,  median  (Q2),  and  the  upper  quartile,  and  the  kernel  density  estimation. 


In  Fig.  16,  relative  differences  of  irradiation  incident  on  tilted 
planes  reach  values  of  10%  for  fixed  systems  and  -10%  for  one- 
axis  and  two-axis  with  higher  dispersion  in  these  last  cases.  KED 
shows  higher  values  than  CM  SAF  especially  in  the  north  area  for 
fixed  systems,  and  to  a  lesser  extent,  for  one  and  two  axis  around 
the  Pyrenees  area.  The  RMSD  and  MAD  values  are  very  similar  for 
all  trackers  and  higher  than  those  corresponding  to  the  irradiation 
on  the  horizontal  plane  (Table  5). 

One  possible  explanation  for  the  positive  values  of  relative 
differences  existing  in  the  area  of  Pyrenees  would  come  from  the 
influence  of  the  terrain  elevation  on  satellite  methods  [44].  The 
solar  irradiation  dependence  with  altitude  is  not  well  described  in 
the  satellite  retrieving  methods  yet.  In  a  mountain  area  each  pixel 
cover  an  area  of  very  varying  altitude  and  therefore  the  irradia¬ 
tion  estimations  have  more  uncertainty  than  in  flat  terrains. 
Besides,  satellite  methods  have  difficulties  in  distinguishing  snow 
and  ice  from  clouds  typically  at  mountain  areas.  An  underestima¬ 
tion  of  the  amount  of  water  vapor  and/or  aerosols  produces  an 
overestimation  in  the  values  of  CM  SAF.  Also  the  uncertainties 
associated  with  estimating  the  clouds  cover  plays  an  important 
role  in  divergence  observed  [14].  A  special  version  of  CM  SAF  for 
mountainous  areas  is  being  developed  using  methods  proposed 
by  and  maybe,  after  this  update,  relative  differences  at  those  areas 
will  decrease  [45]. 


5.  Conclusions 

An  analysis  and  comparison  of  daily  and  yearly  solar  irradia¬ 
tion  from  the  satellite  CM  SAF  database  and  a  set  of  301  stations 
from  the  Spanish  SIAR  network  is  performed  using  data  of  2010 
and  2011.  This  analysis  is  completed  with  the  comparison  of  the 
estimations  of  effective  irradiation  incident  on  three  different 
tilted  planes  (fixed,  two  axis  tracking,  north-south  horizontal 
axis)  using  irradiation  from  these  two  data  sources.  Finally,  a 
new  map  of  yearly  values  of  irradiation  both  on  the  horizontal 
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plane  and  on  inclined  planes  is  produced  mixing  both  sources 
with  geostatistical  techniques  (kriging  with  external  drift,  KED). 

The  comparison  between  the  irradiation  values  from  SIAR,  CM 
SAF  and  KED  is  performed  in  the  context  of  the  SIAR  pyran- 
ometers  tolerance  (5%).  The  difference  of  global  irradiation  from 
SIAR  and  CM  SAF  at  71%  of  the  locations  are  inside  the  range  of 
this  pyranometer  uncertainty.  Outside  this  5%  band,  96.5%  of  the 
SIAR  stations  provide  lower  global  irradiation  values  than  CM  SAF. 
The  relative  difference  increases  when  a  tracking  system  is 
considered:  both  the  standard  deviation  of  the  irradiation  values 
{(Jsiar  and  (Tcmsaf )  and  the  statistics  of  the  differences  (RMSD  and 
MAD)  increase  with  the  use  of  tracking  systems. 

The  Mean  Absolute  Difference  (MAD)  between  CM  SAF  and 
SIAR  is  approximately  4%  for  the  irradiation  on  the  horizontal 
plane  and  is  comprised  between  5%  and  6%  for  the  irradiation 
incident  on  the  inclined  planes. 

The  use  of  kriging  with  external  drift  reduces  the  difference 
with  SIAR  and  CM  SAF,  both  for  the  horizontal  plane  and  for 
inclined  planes.  The  MAD  between  KED  and  SIAR,  and  KED  and 
CM  SAF  is  approximately  3%  for  the  irradiation  on  the  horizontal 
plane  and  is  comprised  between  3%  and  4%  for  the  irradiation 
incident  on  the  inclined  planes. 

6.  Software 

The  methods  described  in  this  paper  have  been  implemented 
using  the  free  software  environment  R  [46]  and  several  contrib¬ 
uted  packages,  namely:  gstat  [40]  and  sp  [41]  for  the  geosta¬ 
tistical  analysis;  solaR  [47]  for  the  solar  geometry,  irradiation 
and  PV  energy  calculations;  raster  [48]  for  spatial  data  manip¬ 
ulation  and  analysis,  and  rastervis  [49]  for  spatial  data  visua¬ 
lization  methods. 

The  source  code  is  available  at  https://github.com/oscarperpi 
nan/CMSAF-SIAR. 
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